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Abstract 

We introduce the nested canalyzing depth of a function, which measures the extent 
to which it retains a nested canalyzing structure. We characterize the structure of 
functions with a given depth and compute the expected activities and sensitivities of 
the variables. This analysis quantifies how canalyzation leads to higher stability in 
Boolean networks. It generalizes the notion of nested canalyzing functions (NCFs), 
which are precisely the functions with maximum depth. NCFs have been proposed as 
gene regulatory network models, but their structure is frequently too restrictive and 
they are extremely sparse. We find that functions become decreasingly sensitive to 
input perturbations as the canalyzing depth increases, but exhibit rapidly diminishing 
returns in stability. Additionally, we show that as depth increases, the dynamics 
of networks using these functions quickly approach the critical regime, suggesting 
that real networks exhibit some degree of canalyzing depth, and that NCFs are not 
significantly better than functions of sufficient depth for many applications of the 
modeling and reverse engineering of biological networks. 

1 Introduction 

A large influx of biological data on the cellular level has necessitated the development 
of innovative techniques for modeling the underlying networks that regulate cell activi- 
ties. Several discrete approaches have been proposed, such as Boolean networks [8], logical 
models [T7], and Petri nets [5]. In particular, Boolean networks have emerged as popular 
models for gene regulatory networks [TJ [T2]. However, not all Boolean functions accu- 
rately reflect the behavior of biological systems, and it is imperative to recognize classes 
of functions with biologically relevant properties. One such notable class is the canalyzing 
functions, introduced by KaufTman [8j, whose behavior mirrors biological properties de- 
scribed by Waddington |20j. The dynamics of Boolean networks constructed using these 
functions are of considerable interest when determining their modeling potential. Random 
Boolean networks constructed using such functions have been shown to be more stable 
than networks using general Boolean functions, in the sense that they are insensitive to 
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small perturbations [10J. Karlssona and Hornquist [7] explore the relationship between 
the proportion of canalyzing functions and network dynamics. In [10J, the authors further 
expand the canalyzation concept and introduce the class of nested canalyzing functions 
(NCFs). In [11], networks of NCFs are shown to exhibit stable dynamics. Also, Nikola- 
jewa, et al. [T3] divide NCFs into equivalence classes based on their representation and 
show how the network dynamics are influenced by choice of equivalence class. Nested can- 
alyzing functions have a very restrictive structure and become increasingly sparse as the 
number of input variables increases [Sj. Also, it is possible that not all variables exhibit 
canalzying behavior. Hence, a relaxation of the nested canalyzing structure is necessary. 

In this article, we further explore canalyzation by analyzing functions that retain a 
partially nested canalyzing structure. We quantify the degree to which a function exhibits 
this canalyzing structure by a quantity we call the nested canalyzing depth. Functions of 
depth d generalize the nested canalyzing functions, because NCFs are the special case of 
when d = k, where k is the number of Boolean variables. In Section |3j we demonstrate 
notable properties of these partially nested canalyzing functions, and show that their rep- 
resentation is unique. This leads to a theorem about the structure of functions of depth 
d, which generalizes a result in [B] about NCFs. In Section |4j we compute the expected 
activities and sensitivities of functions given their canalyzing depths, which are extensions 
of results of Shmulevich and Kauffman [18] about activities and sensitivities of canalyzing 
functions. We prove that as canalyzing depth increases, functions become less sensitive 
to perturbations in the input; however, the marginal benefit incurred by adding further 
canalyzing variables sharply decreases. As a result, functions of larger depth provide an 
improvement in sensitivity over general canalyzing functions, but imposing a fully nested 
canalyzing structure provides little benefit over functions of sufficient canalyzing depth. 
Finally, in Section [5j we use Derrida plots to show that dynamics of networks constructed 
using more structured functions rapidly approach the well-known critical regime, whereas 
networks with functions of relatively few nested canalyzing variables remain in the chaotic 
phase. This is in contrast to the findings of Kauffman et al. [11] . but in agreement with re- 
cent work of Peixoto [16J, and it further supports the biological utility of certain canalyzing 
functions. 

2 Nested Canalyzing Depth 

A Boolean function f(x) = f(x\, . . . , x^) is canalyzing if it has a variable Xi for which 
some input x^ = ai implies f(x) = hi for some hi G {0, 1}. In this canalyzing 
variable, the input a« is its canalyzing value, and the output value hi when the 
corresponding canalyzed value. Note that if / is constant, then every variable is trivially 
canalyzing. 

If a canalyzing variable Xi does not receive its canalyzing input a^, then the output is 
some function gi{xi), where Xi = [x,\, . . . , Xi-\,Xi+\, . . . , x^). If gt is constant, Xi is called a 
terminal canalyzing variable of /. Note that for each i ^ j, Xj is then trivially canalyzing 
in gi- 
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If <?j is not constant, we ask whether it too is canalyzing. If so, there is a canalyzing 
variable Xj with canalyzing input aj, and when Xj 7^ aj, the output of / is a function 
gij(xij), which may or may not be canalyzing. Here, Xij denotes x with both Xi and Xj 
omitted. Eventually, this process will terminate when the function g is either constant or 
no longer canalyzing. 

Definition 1. Let f(x\, . . . ,Xk) be a Boolean function. Suppose that for a permutation 
a G Sk, an integer d > and a Boolean function g(x a ^+i), ■ ■ ■ > x a(k))> 
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where either x a ^) is a terminal canalyzing variable (and hence g is constant), or g is non- 
constant and none of the variables x a ^+i),- ■ ■ , £ CT (fc) are canalyzing in g. Then f is said 
to be a partially nested canalyzing function. The integer d is called the active nested 
canalyzing depth of f , and the full nested canalyzing depth of f is d if g is non-constant, 
and k otherwise. The sequence x a n\, . . . , x a ^) is called a canalyzing sequence for f . 

If we speak of simply the "canalyzing depth" or "depth" of a function, we are referring 
to the full nested canalyzing depth. In the next section, we will show that the depth is well- 
defined, i.e., that it does not depend on the choice of a G Sk- The class of nested canalyzing 
functions (NCFs) [6J [10] are precisely those with active depth k. A constant function (all 
2 k entries in the truth table are the same) is not an NCF by the classical definition, but 
changing a single value in the truth table suddenly makes it nested canalyzing. In our 
set-up, both of these functions have full depth k. Constant functions have active depth 
0. For completeness, we will say that a non- canalyzing function has active and full nested 
canalyzing depth 0. 

Canalyzing and nested canalyzing functions have been used in gene network models 
because they possess biologically relevant features [20J. For example, in a gene regulatory 
network, a collection of k genes that affect the expression level of a particular gene can 
be modeled with a /c-variable Boolean function. While it is believed that some of these 
relationships are canalyzing (e.g., if A is expressed, then B is not expressed, regardless of 
the states of the other genes), it is unreasonable to expect that all relevant genes will act in 
a nested canalyzing manner. For instance, transcription factors in gene regulatory networks 
likely display canalyzing behavior, while other proteins do not. Thus, situations will arise in 
which nested canalyzing functions do not fully capture the dynamics of biological systems. 
For example, in the gene regulatory network of the cell cycle in Saccharomyces cerevisiae, 
genes Cdcl4 and Cdc20 are the only canalyzing inputs for the Swi5 gene. The remaining 
function is constant, which does not fit Kauffman's definition of nested canalyzing, nor 
can this gene's behavior be modeled with an NCF [T2]. Thus, when reverse engineering 
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a biological network with partial data, the rigid NCF structure is restrictive and likely 
inappropriate to model the behavior of the system. Also, the number of NCFs becomes 
rapidly sparse in the set of Boolean functions as k increases. For instance, the proportion 
of NCFs in 6 variables is on the order of 10 -15 [6]. Because of this sparsity, it is unlikely 
that a nested canalyzing function fits a given data set. We will show why functions with 
less than full canalyzing depth exhibit nearly identical key features as NCFs with regards 
to activities, sensitivities, and stability, promoting their potential use in biological models. 

Example 1. Consider the function defined as 

f(xi,x 2 ,x 3 ,x 4: ) = x 2 A (-i£i(A(x 3 XOR x 4 ))) , 

which is similar to Example 2.4 in [6]. There is a canalyzing sequence X2,Xi, with a 2 = 
b 2 = ; and ai = 1 and hi = 0. The remaining function g{x3 1 x^) = X3 XOR X4 is not 
canalyzing in either variable. Therefore, f is a PNCF of active and full depth 2. 

3 Properties of Partially Nested Canalyzing Functions 

Proposition 1. Let f(x) be a k-variable Boolean function. Then 

(i) If Xi = a, implies f(x) = hi, and x i: 7^ a { leaves gi(xi), then at least half of the truth 
table values of f must be b{. 

(ii) If exactly half the truth table values of f are h i; then either Xi is terminally canalyzing, 
or f is a non- canalyzing function. 

Proof. The statement in (i) follows because Xi = for exactly half of the input values in 
the truth table. The corresponding output value must be hi for at least these inputs. To 
show (ii), suppose that / is canalyzing, and = implies f(x) = hi. By (i), 7^ 
implies f(x) = ->bi. Therefore, g(xi) is constant and Xi is terminally canalyzing. □ 

The canalyzing depth of a function can be computed in a divide-and-conquer manner 
described in Algorithm [TJ The algorithm scans the truth table for a canalyzing variable, 
and upon finding one, removes the columns for which the canalyzing variable takes the 
canalyzing input value. This is repeated until no more canalyzing variables are present or 
a constant function remains. Proposition [T] and the structure of the truth table imply that 
if there is a tie for b, then there is a terminally canalzying variable or there are no canalyzing 
variables. Therefore, it is not necessary to test both b and ->b as possible canalyzed values. 
In the execution of the algorithm, we set a flag whenever a tie for canalyzed value arises. 
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Algorithm 1. 

1. Setd = 0. For i = 1. . . k - 1 : 

(a) Set 6 = 0, flag = 0. Let I be the number of ones in the truth table. 

• If £ == 2 k ~ l+1 return k. // Constant function remains 

• If £ == 2 k ~ % , set flag = 1. // Tie in output value 

• If£>2 k -\ b=l. 

(b) For remaining k — i + 1 variables in truth table: 

i. Let x be the number of input ones and y the number of input zeros that give 
output b. 

ii. * If x == 2 k ~ % , the current variable is canalyzing with input 1 and output 

b. Remove canalyzing rows and current variable from truth table and 
break out of loop. 

• If y == 2 k ~ % , the current variable is canalyzing with input and output 
b. Remove canalyzing rows and current variable from truth table and 
break out of loop. 

(c) If no variables were found to be canalyzing, return d; else d++. 

(d) If flag == 1, return k. // Constant function remains 

2. Return k. 



Note that it takes exponential time simply to view the entire truth table of /; however, 
the algorithm is linear in the size of the table. Indeed, the i th step of Algorithm [l] takes 
(k — i) • 2 k ~ l steps, and so the running time is 

J2(k ~ ^)2 k - l <k-2 k (l+ ] - + - + ..\=0{k-2 k ). 
i=i ^ ' 

We can use Algorithm [T] to show that canalyzing depth is well-defined, meaning that it 
does not depend on the choice of canalyzing sequence. First, in Theorem [T] we present the 
functional form of a PNCF of canalyzing depth d. Lemma 2.6 in [6] is the special case of 
Theorem [T] when / is nested canalyzing, i.e., when d = k. 

Theorem 1. Let yi = x a ^) + ai + bi, 1 < i < d and let 

f(xi, ...,x k ) = yiOi(2/a02(- • • (Vd<>d9(x <7{d+1) , . . . , x a{k) )) ...)), (2) 

J V if b t = l 
1 A if h = ' 

(i) None of the variables x ff u+i), ■ ■ ■ ,Xo-(k) are canalyzing in g, or 

(ii) g is a constant function. 



where 



0i 



di, bi e {0, 1} for 1 < i < d, and 
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Then f has canalyzing depth d, with canalyzing sequence x a m, . . . ,x a u). These variables 
have canalyzing values a%, . . . , ad and canalyzed values b±, . . . ,bd- Furthermore, any func- 
tion of canalyzing depth d can be represented in this form. 

Proposition [T] and our previous observations indicate that in case of a tie in potential 
canalyzed values, we cannot make a "wrong" choice for b in the execution of Algorithm [T] 
Hence, to show that the depth is unique, it suffices to show that if there are multiple 
canalyzing variables at a given iteration, the depth does not depend on our choice of 
canalyzing variable. 

Proposition 2. The nested canalyzing depth computed using Algorithm^ yields a unique 
answer. 

Proof. Suppose that at a given iteration there are m variables left in the truth table, and 
two of them are canalyzing. Note of the 2 m values in the truth table, there are 2 m_1 
canalyzing inputs for one variable and 2 m_1 for the other, 2 m ~ 2 of which are canalyzing 
for both. Therefore, by Part (i) of Proposition [TJ the two canalyzed output values must 
be the same. Also, regardless of which variable enters at the current iteration, the other 
variable will have 2 m_2 canalyzing input values remaining at the next iteration, and will 
still be a canalyzing variable at the next iteration with the same canalyzed value. Now, 
notice in Theorem [T] above that if 0« = Oi+i for i = 1, . . . , k — 1, interchanging the variable 
order for x a (^ and Xo-^+i) does not change the function. Since <\ is determined by the 
canalyzed output the function is the same regardless of which variable is selected 
at the current iteration. An analogous argument holds when more than two canalyzing 
variables are present. □ 

It is now easy to see that the nested canalyzing structure introduced in Equation [T] is 
well-defined since the remaining function g is unique. 



4 Activities and Sensitivities 

In this section we compute the expected activities and sensitivities of functions based on 
their canalyzing depth, and in the next section we will tie these results to the stability 
of Boolean networks based on the canalyzing depth of the individual functions. Let x G 
{0, l} k , and write x- 7 '* = (x\, . . . , Xj-i,i, Xj + i, . . . , Xk) and let © be the XOR function. The 
partial derivative of f{x\, . . . , Xk) with respect to Xj is 

^ = /(x-°) © /(x- 1 ). 

The activity (or influence) of a variable Xj in / is 



. 1 \^ #/(x) /ox 

^^ = 9* 2^ ( 3 ) 



•' l ^ dXj 

<e{o,i} fe 
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and the sensitivity of / is denned by 

k 

S '(x) = ^ X [/(x® e ^/(x)], 
1=1 

where is the i th unit vector and \ is an indicator function. The activity a? quantifies how 
often toggling the j th bit of x toggles the output of /, and the sensitivity s* (x) measures 
the number of ways that toggling a bit of x toggles the output of /. The average sensitivity 
of / is the expected value of s* (x) taken uniformly over all x G {0, l} k , i.e., 

k 

st = E[ S f{x)]=Y,oi. (4) 

i=l 

In [IB], Shmulevich and Kauffman show that a random unbiased Boolean function in 
k variables has average sensitivity |. Also, they prove that for an unbiased canalyzing 
function (i.e., depth at least 1) with canalyzing variable X\, the expected activities of 
(xi, . . . , Xk) are given by 

and hence the average sensitivity is = ^j-k The following theorem extends this to 
functions of arbitrary canalyzing depth. 

Theorem 2. Let f be a Boolean function in k variables with nested canalyzing depth at 
least d. Renumbering the variables if necessary, assume that x\,...,Xd is a canalyzing 
sequence. Then, if we assume a uniform distribution on the function inputs, the expected 
activities of the variables (x\, . . . , Xf.) are given by 




Furthermore, the expected sensitivity of f is 

i=l 

Proof. Since we are assuming a uniform distribution on the function inputs, for any variable 
arj, 1 < j 1 < k, we can think of the activity of Xj as the probability that changing the input 
to the j th entry changes the function output. That is, 

af M) = ^ E ^ = W(x©e,)^/(x)). 
xe{o,i} fe 3 
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Now, if Xj is a canalyzing variable, we know by Equation [TJthat if at least one of x±, . . . , Xj_i 
gets its canalyzing input, the input to Xj cannot affect the function output and this prob- 
ability is 0. Hence, we have 

af = P(/(x© ei )^/(x)) 

= P(/(x © ej) ^ /(x)|xi 7^ ai, . . .,xj-i ^ a j ^ 1 )P(xi ^ a h . . . ,Xj-i ^ 

Since each canalyzing variable receives its canalyzing input with probability |, 

P(xx ^ a u . . . , Xj-i ^ a,_i) = 
Also, since / is a random, unbiased function, 

P(/(x©e i ) ^ /(x)|a;i ^ oi,...,^! ^ a^i) = -. 

Therefore, 

/ _ 1 1 1 
a J ~ 2 ' 2^1 ~ 

Alternatively, if Xj is a non-canalyzing variable, the input to Xj is only relevant when 
none of the canalyzing variables X\,. . . ,x d get their canalyzing inputs. Using a similar 
argument as above, we see that 

aJ = P(/(x©e,)^/(x)) 

= P(/(xffie i ) ^/(x)|xi ^ a l ,...,x d ^ a d )P{x x ^ a 1 , . . . , x d ^ a d ) 
1 1 1 



2 2 d 2 d+1 ' 

Equation [7] now follows from Equation |4j □ 

Note that this theorem may also be proven via induction on d, with Equation [5] as a 
base case, following an argument similar to that in [T8] . 

By Theorem [2j the average sensitivity of a function decreases as the depth increases. 
However, the differences in sensitivity become increasingly smaller, and are precisely 

E[sf«] - E[sf*»] = 1- - + — -! + — - 
k-d-l 

= — ^d+2 — — ^ ' when k — d > 1 . 

Observe that this quantity rapidly goes to zero, and so each subsequent canalyzing variable 
has a much smaller impact on the sensitivity. Thus, the difference in sensitivity between 
fully nested canalyzing functions and partially nested canalyzing functions of sufficient 
depth is very slight. For example, Table |T] gives the expected sensitivities for PNCFs with 
k = 6 and k = 12 input variables, respectively. To compare the two, we normalize the 
depth d by the number of variables k. 
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Table I: The expected sensitivity E [s^] for PNCFs in k variables of depth d. 



d/k 





1/6 


1/3 


1/2 


2/3 


5/6 


1 


k = 6 


3.0000 


1.7500 


1.2500 


1.0625 


1.0000 


0.9844 


0.9844 


k — 12 


6.0000 


2.0000 


1.1875 


1.0313 


1.0039 


1.0000 


0.9998 



5 Stability and Criticality vs. Canalyzing Depth 

Boolean networks created using classes of functions with a lower sensitivity have been 
shown to be more dynamically ordered than those with a higher sensitivity [IB]. This 
stability is an important feature of biologically relevant functions, and so it is essential to 
determining the utility of such functions as biological models. In order to quantify the 
extent to which functions with larger depth (and hence smaller sensitivity) result in more 
dynamically stable Boolean networks, we constructed random Boolean networks composed 
of PNCFs of varying depth. We used the annealed approximation mean-field theory due 
to [3] and Derrida curves to display the results. The curves are defined as follows. Let 
x 1 (i) and x 2 (t) be two states in a random Boolean network, and define p(t) to be the 
normalized Hamming distance, i.e., p{t) = - ■ ||x 1 (t) — x 2 (i)||i, where || • ||i is the standard 
i 1 metric. The Derrida curve is a plot of p(t + 1) versus p(t) averaged uniformly over 
different states and networks. If the curve for small values of p(t) lies below the line 
y = x, then small perturbations are likely to die out, and the network is said to be in 
the frozen phase. The phase spaces of frozen networks consist of many fixed points and 
small attractor cycles. If the curve lies above the line y = x, then small perturbations 
generally propagate throughout the network, and the network is said to be in the chaotic 
phase, characterized by long attractor cycles. The boundary threshold between these two 
is the well-known critical phase [I] . It has been recently suggested [21 HH 151 EES] that many 
biological networks tend to lie in the critical phase, as these systems must be stable enough 
to endure changes to their environment, yet flexible enough to adapt when necessary. 

We constructed ensembles of randomly wired networks with n = 100 nodes, each with a 
randomly chosen Boolean function with k = 12 variables. We chose the individual functions 
by sampling uniformly across the class of PNCFs of depth at least d, for d = 0, 2, 4, . . . , 12. 
We will refer to such a network as a depth-d network. To sample uniformly across all PNCFs 
of depth at least d, we used a random number generator to select d nested canalyzing 
variables, and a permutation a of these variables. We then used a random bit generator to 
select the canalyzing values ai, . . . , a<j and canalyzed values bi, . . . bj. We had a potential 
bias in our function selection arising from the fact that if 0i — 0i+i (or equivalently, 
bi = b i+ i, as bi determines <>«), then interchanging a{i) with a{i + 1) does not change the 
function. To eliminate this bias, we only allowed functions where o~(i — 1) < a(i) whenever 
0,;_i = 0i for i — 2, . . . ,d — 1. Finally, we used a random bit generator to determine the 
function in the remaining k — d variables. Our sampling method for creating the random 
networks is similar to [IB]. For each d, we created 25 random Boolean networks using 
functions of said depth and sampled from each network. Since p(t + 1) is determined 
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Derrida Curve with n=100, k=12 
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Figure 1: Derrida curves for random Boolean networks with n = 100 nodes and k = 12 
inputs per function. 



experimentally, we computed it as the sample mean, sampled over the depth-c? random 
networks for each depth. We also constructed Derrida curves using the sampling method 
described in [9] , which generated nearly identical results. The resulting Derrida curves are 
shown in Figure [T} 

The Derrida curves corresponding to networks constructed using functions of larger 
depth show more orderly dynamics than those of smaller depth. This reaffirms the idea 
in [IB] that sensitivity of a function is an indicator of the dynamical stability of networks 
constructed with these functions. The curves move closer together as the depth increases. 
For example, the depth-2 networks are much more stable than the depth-0 networks, and 
networks with functions of depth at least 4 are even more stable; however, the marginal 
benefit of stability as depth increases drops off sharply - the Derrida curves are nearly 
identical for networks with functions of depth 4, 6, 8, 10, and 12. This matches our 
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theoretical results of Theorem [2] on expected activities and sensitivities and illustrates how 
the earlier canalyzing variables have a much greater influence. It also suggests that there 
is little benefit in imposing the full nested canalyzing structure in network models, as 
functions with large enough canalyzing depth exhibit very similar stability results without 
the rigidity of being fully nested canalyzing. Additionally, for small values of p(t), the 
curves quickly approach the line y = x from above, indicating that these networks rapidly 
move from the chaotic phase toward the critical phase. This is contrary to the claim of pU] 
that networks comprised of canalyzing functions are always in the frozen phase, but it is 
in alignment with recent findings of [16] which also refute Kauffuman's claim, accrediting 
his results to his choice in parametrization. 

6 Concluding Remarks 

Canalyzing and nested canalyzing functions have been proposed as gene network models 
because they exhibit biologically relevant properties. While it is reasonable to expect some 
Boolean models to have functions with some degree of nested canalyzation, fitting biological 
data to fully nesting canalyzing functions can be at times artificial, and at other times 
simply incorrect. Our analysis of the depth that a function retains a canalyzing structure 
elucidates the role of canalyzation in the dynamics of networks over these functions. Our 
results on the structure of PNCFs generalize known results on NCFs, and our results 
on the activities and sensitivities of variables in functions of a given depth generalize 
similar theorems of simple canalyzing functions. Moreover, we saw that in random Boolean 
networks, the stability increases with canalyzing depth. However, the marginal gain in 
stability drops off quickly, in that the stability of our networks with functions of depth 
at least d — | were nearly identical to those with full depth d = k. Additionally, just a 
few degrees of canalyzation are necessary to drop the network into the critical regime, in 
which many real networks are believed to exist. Together, this suggests that using NCFs 
in biological models for stability reasons is not only at times rather contrived, but simply 
unnecessary. 
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